#ifndef _inner_product_hpp_
#define _inner_product_hpp_

#include <numeric>

#include "vector_traits.hpp"

namespace core {

    namespace math {

        namespace detail {

            template<class V, class F1, class F2, std::size_t n>
                struct inner_product_unroll_loop
                {
                    typedef typename field<V>::type  field_type;
                    static field_type eval(const V& v1, const V& v2, const field_type& init_value, F1 binary_op1, F2 binary_op2)
                    {
                        return binary_op1(inner_product_unroll_loop<V,F1,F2,n-1>::eval(v1,v2,init_value,binary_op1,binary_op2),binary_op2(v1[n],v2[n]));
                    }
                };

            // break condition
            template<class V, class F1, class F2>
                struct inner_product_unroll_loop<V,F1,F2,0>
                {
                    typedef typename field<V>::type  field_type;
                    static field_type eval(const V& v1, const V& v2, const field_type& init_value, F1 binary_op1, F2 binary_op2)
                    {
                        return binary_op1(init_value,binary_op2(v1[0],v2[0]));
                    }
                };

            struct inner_product_static_impl
            {
                template<class V, class F1, class F2>
                    static typename field<V>::type eval(const V& v1, const V& v2, const typename field<V>::type& init_value, F1 binary_op1, F2 binary_op2)
                    {
                        return inner_product_unroll_loop<V,F1,F2,
                               dimension<V>::value-1 >::eval(v1,v2,init_value,binary_op1,binary_op2);
                    };
            };

            struct inner_product_dynamic_impl
            {
                template<class V, class F1, class F2>
                    static typename field<V>::type eval(const V& v1, const V& v2, const typename field<V>::type& init_value, F1 binary_op1, F2 binary_op2)
                    {
                        std::size_t dim =   dimension<V>::get(v1);
                        //check that all both vectors have the same dimension
                        assert( dim == dimension<V>::get(v2));
                        using std::inner_product;
                        return inner_product(v1.begin(),v1.end(),v2.begin(),typename field<V>::type(0),binary_op1,binary_op2);
                    };
            };

            struct inner_product_scalar_impl
            {
                template<class V, class F1, class F2>
                    static V eval(const V& v1, const V& v2, const V& init_value, F1 binary_op1, F2 binary_op2)
                    {
                        return binary_op1(init_value,binary_op2(v1,v2));
                    };
            };

            template<class V, class F1, class F2>
                struct choose_inner_product_impl
                {
                    //check if this is a vector with static dimension
                    typedef typename boost::mpl::eval_if<is_scalar<V>,
                                                boost::mpl::identity<inner_product_scalar_impl>,
                                                boost::mpl::eval_if_c<(dimension<V>::value>0),
                                                                    boost::mpl::identity<inner_product_static_impl>,
                                                                    boost::mpl::identity<inner_product_dynamic_impl>
                                                                    >
                                                >::type  impl_type;
                    typedef typename field<V>::type scalar_type;
                    static scalar_type eval(const V& v1, const V& v2, const scalar_type& init_value, F1 binary_op1, F2 binary_op2)
                    {
                        return impl_type::template eval<V,F1,F2>(v1,v2,init_value,binary_op1,binary_op2);
                    }
                };

        } // namespace detail

        template<class Vector, class BinaryFunctor1, class BinaryFunctor2>
            typename field<Vector>::type inner_product(const Vector& v1, const Vector& v2, const typename field<Vector>::type& init_value,
                    BinaryFunctor1& binary_op1, BinaryFunctor2& binary_op2)
            {
                typedef detail::choose_inner_product_impl<Vector,BinaryFunctor1&,BinaryFunctor2&> impl_type;
                return impl_type::eval(v1,v2,init_value,binary_op1,binary_op2);
            };

        template<class Vector, class BinaryFunctor1, class BinaryFunctor2>
            typename field<Vector>::type inner_product(const Vector& v1, const Vector& v2, const typename field<Vector>::type& init_value,
                    const BinaryFunctor1& i_binary_op1, const BinaryFunctor2& i_binary_op2)
            {
                BinaryFunctor1  binary_op1(i_binary_op1);
                BinaryFunctor2  binary_op2(i_binary_op2);
                typedef detail::choose_inner_product_impl<Vector,BinaryFunctor1&,BinaryFunctor2&> impl_type;
                return impl_type::eval(v1,v2,init_value,binary_op1,binary_op2);
            };

    } // namespace math

} // namespace core

#endif // _inner_product_hpp_
